Decaying of Phase Synchronization - A Physiological Tool 
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We describe the effects of the asymmetry of cycles and 
non-stationarity in time series on the phase synchronization 
method. We develop a modified method that overcomes these 
effects and apply this method to study parkinsonian tremor. 
Our results indicate that there is synchronization between 
two different hands and provide information about the time 
delay separating their dynamics. These findings suggest that 
this method may be useful for detecting and quantifying weak 
synchronization between two non-stationary signals. 



INTRODUCTION 

Analysis of the synchronization between subsystems of 
a complex dynamical system is important for character- 
izing the system. Much effort has recently been focused 
on detecting and quantifying synchronization within non- 
stationary noisy systems (see e.g., [1,2]). In such systems, 
the traditional cross-correlation technique may not assure 
appropriate detection in all cases of interdependency [2] . 
The phase synchronization method has been found useful 
in identifying synchronization in several systems. It has 
been used to detect brain activity associated with parkin- 
sonian tremor [3], to describe cardiorespiratory systems 
[4,5], and to find relations between temperature and pre- 
cipitation in different regions [6]. 

Here we further develop the synchronization decay 
method and apply it to clarify some questions related to 
parkinsonian tremor. Using the accelerometric and elec- 
tromyographic (EMG) records of the hands during "rest 
tremor" , we show that the two hands' tremor movements 
are synchronized and that the coupling mechanism has 
a broad range of time response scales starting at scales 
that are below the detection of the measuring apparatus 
(the sampling rate is 400Hz) and reaching 0.5 seconds 
and above. 



EFFECTS OF ASYMMETRY AND 
NON-STATIONARITY 

Direct application of the phase synchronization 
method [2] may encounter difficulties regarding consis- 
tent asymmetry of the density of points in each cycle, as 
discussed below. The phase synchronization method uses 
the concept of "analytic signal" (AS) [8] which is useful 



for detecting the phase in periodic and semi-periodic sig- 
nals [2]. By computing the AS S A (t) = S(t) + iS h (t) 
(where S is the original signal, and Sh is its Hilbert trans- 
form, see formulation in the method section) of two si- 
multaneous records S a (t) and Sb(t), and evaluating the 
phase differences series ip{t) = <j> a (t) — 4>b{t) (where the 
phase 4>a,b is the angle of the AS in complex space) one 
can generally estimate how strongly the signals are cou- 
pled, observing the variation of ip. However, if the ana- 
lytic signals are not folded equally over the cycle , their 
phases tend to stay longer in the more folded areas (see 
Fig. la, for an example of asymmetric folding). In this 
case, synchronized dynamics, where the completion of a 
cycle in S a depends on the completion of a cycle in Sb, 
would be hard to detect. This difficulty arises because S a 
may complete a very folded area (and remain in a very 
narrow region in phase space) while Sb completes a much 
less folded area (and therefore swaps a wide angle in the 
complex space). Thus, the phase differences series ip [1] 
will not have a very dominant phase lag, and the synchro- 
nization might appear to be insignificant. The geometry 
of the AS might have, therefore, a significant impact on 
the ability to detect synchronization even when the dy- 
namics of S a and Sb are coupled. 

Non-stationarity of S a and Sb may lead to another dif- 
ficulty. It is accepted that the global average of the AS 
must be removed before extracting the phase series, oth- 
erwise, the AS will appear to be "folded asymmetrically" 
even if it is a perfect circle (because the circle is not 
centered at the origin). If the center of the AS ring is 
non-stationary, this would be the case even after one re- 
moves the global average. A partial solution will be to 
identify small regimes, where the center of the ring does 
not wonder too far, remove the average (detrending), and 
calculate the synchronization index in each regime. Fi- 
nally we average the synchronization index [3] over all 
regions. This compromise may still imply quite small 
synchronization indices when the signal is highly non- 
stationary. Non-stationarity will therefore also have high 
influence on the synchronization index. 

The original work of Huygens [9] discussed the syn- 
chronization of pendulums due to interactions through 
the embedded media. In this case, a typical response 
time exists. The dynamics of one oscillator will gener- 
ally follow the other oscillator in a time delayed manner. 
One, then, should expect that choosing the right time 
delay r will improve the Synchronization index of the 
signals S a (t) and Sb(t-r). 

We therefore describe a procedure (also used partially 
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in [6]) with a related new index, which detects interde- 
pcndcncics between processes, taking into account asym- 
metry and non-stationarity, and providing a way to es- 
timate the typical time delay [7]. We will assume that 
interdependencies can be observed from phase differences 
alone, and that the process of mutual modification does 
not repeat itself exactly but has a stochastic component. 



THE METHOD 

In order to evaluate the significance of the synchroniza- 
tion indices measured, we compare it with the "natural" 
values of the index obtained when the signals arc not 
coupled. We therefore use the concept of "the synchro- 
nization decay" [6] . This concept assumes that since the 
coupling is local in time, finding the synchronization of 
the first time series with the future (past) of the second 
time series will give a synchronization value that is re- 
lated to the geometric structure of the signal, and not to 
the influence of the signals on each other. The decay of 
synchronization is, therefore, just a plot of the "simple 
synchronization" indices (as defined below) as a function 
of the time differences between the records. The algo- 
rithm we use to estimate the decay of synchronization is 
as follows: 

Firstly we choose a reasonable time shift T max within 
which we can assume the local influence of one oscillator 
on the other is negligible. We set the time shift r to be 
—Tmax- If T max = and only one window (see below) 
is taken, the method is equivalent to the standard phase 
synchronization method. 

1. Take the signals Si(t), S2(t — r) when 
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2. Evaluate the analytic signals 



sr A (t) = sr(t) + isr H (t) 



S-(t) + — *S T (q) 
irq 



where the superscript r represents the time shift 
that is taken between the signals, and "*" is a con- 
volution with the integration variable q. Evaluation 
is done using the AS form in the frequency domain 



/ dte-^(s(t)+ijdq£^) = 

dt dq e-^S{q) (<J(i - q) + ^) - 
(6(q) + ^)*S(q) = (l + sgn(f))S(f) = 
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3. The analytic signals can be represented in the form 
S A {t) = A(t)exp(i(t) T (t)). The phase series <j) T (t) 
can be easily obtained by the trigonometric relation 



(t) = arctan 



S T „(t) 



+ -(l- sgn (S T )). 



4. Calculate the series 

CmO) =mod2ir m4>l{t) - n<j>l(t) 

where m,n are integers (^ 0). The ratio m:n equals 
to the ratio between the dominant frequencies uj 2 ■ 

UJl. 

5. Evaluate the probability density P of ^ mn {t) , by 
choosing a fine division of [0, 27r] to bins [10], 
and evaluating the probability Pi of finding ip^n (*) 
within the i-th bin. [15] 

6. Evaluate the Shannon entropy [11] of the probabil- 
ity histogram 



S T = -J2PilogPi 



Derive the synchronization index p [3] by scaling 
the entropy S T to be in [0, 1] where stands for 
lack of synchronization and 1 stands for full syn- 
chronization, using the relation 

C _ or 
T ^max ^ 

9 = ^ 

^max 

where S max = logiV. 

7. Set r = t + 1, return to step 1. 

The algorithm should be applied on relatively short time 
scales so that the slips of the center of the AS' ring can be 
ignored. Nevertheless the time window should be wide 
enough to contain many cycles (so that edge numerical 
effects can be ignored). The decay of the synchroniza- 
tion index for all time windows should be averaged over 
all windows. A detailed picture of the development of 
coupling in time can be achieved by observing the indi- 
vidual time windows. 

A "synchronized state" is a local state. It occurs when 
the synchronization index vs. time shift r has two sig- 
nificantly falling tails and one or a few maxima (e.g., 
see Figure 2). A well pronounced synchronization occurs 
when the difference between the tails and the maxima is 
more than the standard deviation of the tails. We quan- 
tify the significance by evaluation of 

< Pke middle > — < Pke tails > 
\j < Pketails > - < Phe tails > 2 
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APPLICATION TO TREMOR IN PARKINSON'S 
DISEASE 

Parkinsonian rest tremor acceleromctric records, being 
non-stable, gave rise to an interesting question, where the 
synchronization decay diagrams might give some indica- 
tion for the answer. We analyzed 29 records, from 3 male 
and 2 female patients, with ages between 55 and 69. Each 
record has a duration of 330 seconds, and a sampling rate 
of 400Hz. 

The Parkinsonian rest tremor accelerometric and EMG 
records in different extremities have been found to be 
non coherent [13], despite their concurrent motion. It 
can be shown that the low coherence between two time 
series imposes very low Synchronization Index as well. 
" Synchronization Decay" diagrams can indicate that al- 
though the Synchronization Indices are indeed low, the 
hands' dynamics might be coupled. We find that syn- 
chronization between the hands appears to be maximal 
with a time delay which suggest a non mechanical cou- 
pling (assuming that mechanical coupling time scale is 
very different from neuronal coupling time scale). 

We tested whether coupling between the dynamics of 
the extremities exists. There is evidence in human hand 
tremor [12] and in neuronal activity of MPTP monkey 
models [14] that the first two modes of oscillations (e.g., 
4 and 8 Hz) exhibit different dynamics which can not 
always be described as harmonics. We therefore com- 
puted the synchronization decay between the hands and 
between each of the modes of the left hand and each of 
the modes of the right hand. When asserting phase syn- 
chronization decay between the two full signals, one can 
refer to the phases extracted as "effective phases", and 
to the phases of the individual modes as "phases". We 
attribute interdependence for cases in which the decay 
of synchronization is larger than 1.5 standard deviations 
of the tails. We found that out of 29 records of the en- 
tire signal, 19 were synchronized. The first mode of the 
tremor hand was synchronized with the first mode of the 
relatively healthy hand in 16 of the records, and with 
its second mode in 11 cases. The second mode of the 
tremor hand was synchronized with the first mode of the 
relatively healthy hand in 6 cases, and with its second 
mode in 7 cases. These results suggest that different ex- 
tremities oscillate dependently. A close look at Fig 3 re- 
veals that synchronization of the first mode of the tremor 
hand and any of the modes of the relatively healthy hand 
is in some cases more sharp (above 3 standard devia- 
tions) than the synchronization of the second mode of 
the tremor hand with any of the modes of the relatively 
healthy hand. The appearance of a significant interde- 
pendence between oscillators, when their coherence is in- 
significant [13], is a surprising result of the analysis. The 
other possibility, that coherence can exist without syn- 
chronization, has been recognized before [2]. The reason 
for this new possibility is due to our different (less strict) 
demand for observing inter-dependency. We suggest that 



it is sufficient to demand that synchronization decays as 
time shift grows. 

We assume above that our findings are related to the 
disease state, and that the outcoming mechanical oscilla- 
tion is indeed a reflection of the neuronal activity. How- 
ever, an alternative explanation is that the mechanical 
system itself (the hands and the body) causes the rela- 
tions found and not the pathology. Thus, we must ex- 
clude the possibility that the coupling we tracked can be 
fully explained by the mechanical coupling between the 
hands through the body media. Our suggestion that it 
is due to neuronal activity is supported by three results: 
(a.) The range of response times (between times less than 
2 msec and up to 0.5 sec) that was found is recognized 
as typical neuronal time scale, (b.) In two of the cases, 
there was an order of magnitude difference between the 
response time of second mode (in both hands) and the 
response time between the first mode in one hand and the 
second mode in the other. Thus, even the same person, 
with the same mechanical parameters shows extremely 
different response times (which also supports the claim 
that there are indeed several sources which rule each 
hand, see [12]), and therefore , a mechanical coupling 
can not explain such delays, (c.) We also tested the 
synchronization decay between the EMG signals. The 
EMG records are quite noisy, and in most cases analysis 
gave no sign of the synchronization such as that found 
between the accelerometric records. However there were 
6 (out of 29) records which did show good synchroniza- 
tion decay. Such interdependence is better explained by 
the interdependence of neuronal activity rather than me- 
chanical couplings (because we measure the signal from 
the brain, and not the actual movement). 

DISCUSSION 

A. Oscillations in the synchronization decay plot 

It is worth noting that one can usually identify the arti- 
facts of non-stationarity by looking at a synchronization 
decay plot. For example, in Figure 2b one can see that 
the synchronization index does fall dramatically. The 
falling profile is modulated with an oscillation. This os- 
cillation can be explained with a reasoning similar to that 
we gave for the index falling due to geometry in the dis- 
cussion above. When shifting the records so that folded 
areas fall together most of the time, one obtains "bet- 
ter" synchronization. When folded area of one record 
overlaps the unfolded areas of the other record (and vice 
versa, of course), the index falls. Therefore the index 
falls and grows alternately as we shift the records (in the 
process described above). Oscillation, is therefore, due 
to geometry. 

B. Why two perfect sines are not synchronized? 

Synchronization of two sines (or of two pendulums with 
the same rod length), using our new general restriction 
(that existence of synchronization is indicated by the fact 
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that the cycles of one oscillator are less coherent with the 
future cycles of the other oscillator than they are with its 
present), is very low (because there will be no change in 
the synchronization index between different time shifts). 
Synchronization means therefore a process of local (in 
time) feedbacks. 

A simple illustration for this point is given in Fig. 4. 
A triangular wave is given by: 

2h ( 00 t 

«(*) = - EM-* + (i + 2i)-)- 

sgn(-t + ir)](t - (0.5 + 2i)^)+ 

[agn(-t+(2 + 2i)^)- 

sgn(-t + (1 + 2i)I)](-i + (1.5 + 2i) T -)} , 

where ^ is the amplitude, and r is the period. A decay 

of synchronization index between the two first modes of 
oscillation of a perfect triangular periodic wave does not 
exist because each cycle of the first mode is prefectly syn- 
chronized with each cycle of the second mode, and not 
just with the relatively close cycles (r is not a function 
of time, therefore each cycle is identical in its phase dy- 
namics with previous cycles). Each point on the decay 
of synchronization plot is a local maximum. 

Even if you perturb by adding (u'(t) = u(t) + 677(f) 
when 77(f) is some white noise) or multiplying (h = h(t)) 
the wave by some random noise, you won't get decay of 
the synchronization index, because phases stay the same 
(detection of the phases might differ a bit, but not consis- 
tently enough to reduce significantly the Synchronization 
Index). In Fig. 4a,b,c there is a superposition of several 
noisy triangular waves (r is fixed for each of them, and 
h = h(t)). It's obvious that the Synchronization Index 
is not 1 all over the diagram (because detecting is not 
perfect) but there's no decaying (Fig. 4c). Nevertheless, 
perturbing the frequency [16] (t = r(f)Fig. 4d,e,f) mim- 
ics a dependency between the two modes of oscillation, 
because a local change in the proceeding of one mode is 
seen simultaneously in the other mode. In real systems, 
such "simultaneous" local changes are the heart of a dy- 
namical synchronization process. This simple example 
suggests that harmonies of a signal with changing local 
period synchronize (in the sense of decaying synchroniza- 
tion) in order to build the wave form. Indeed, the decay- 
ing of the synchronization index plot in Fig. 4f shows 
sharp decay in both tails. Simultaneous local (in time) 
changes in two time series (and therefore - the measure 
of synchronization decay) indicate interdependence. 
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FIG. 1. Example of a parkinsonian tremor record's ana- 
lytic signal as a function of time. The analytic signal always 
oscillates in a counter clock wise direction, 
a: the full signal's transform. 

b: analytic signal of only first mode of oscillation. 
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FIG. 2. (a) Synchronization decay between oscillations in 
both hands. Records here are not filtered, (b) Synchroniza- 
tion decay of first mode in both hands. There are clearly 
decaying tails, but the figure is modulated by a strong oscil- 
lation. In both cases one sees that p , the standard synchro- 
nization index, is extremely small, but is significantly larger 
than the background noise. 
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FIG. 3. Distribution of difference between center and tails 
mean values (in standard deviation units) in the synchroniza- 
tion decay diagram testing interdependency between a. First 
mode of healthy hand vs. first mode of tremor hand. b. Sec- 
ond mode of healthy hand vs. first mode of tremor hand. c. 
First mode of healthy hand vs. second mode of tremor hand, 
d. Second mode of healthy hand vs. second mode of tremor 
hand. e. Full Tremor signals of both hands, f. Full EMG 
signals of both hands. 
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FIG. 4. a,b,c : A superposition (with random in time co- 
efficients) of 10 triangular waves each with slightly different 
wave length, the separation into two dominant modes and the 
Decay of the Synchronization Index plot. 
d,e,f : A triangular wave, its separation to modes and the De- 
cay of the Synchronization index plot. On each cycle a new 
wavelength is chosen randomly. 
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